*! by Hung-Jen Wang
*!    Department of Economics, National Taiwan University
*!    wangh@ntu.edu.tw

*! Compute the imaginary erf function (i.e., erfi).
*! erfi(x) = (2/sqrt(_pi)) * int_0^x exp(t^2) dt

*! Update on Aug 17, 2010.
*! Found that the program returns missing value when "upper" is large (> 9.5379).
*! So set a limit to it.

capture program drop myerfi
program define myerfi

syntax , UPper(string) GENerate(string)

* upper: the upper limit of the integral
* generate: the value is to save in the name of the scalar

* ex: myerfi, upper(1.5) generate(myvalue)

preserve
tempfile temp1

tempvar x y Sy

if `upper' >= 9.5379 {
  scalar `generate' =  1.918e+38
}


else {

quie range `x' 0 `upper' 1000 /* I found 1000 work sufficiently good */


quie generate double `y' = exp((`x')^2)


quie integ `y' `x', gen(`Sy')

scalar `generate' = r(integral)*2/sqrt(_pi)

}

restore

end
